Hopping Transport in Hostile Reaction-Diffusion Systems 
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We investigate transport in a disordered reaction-diffusion (RD) model consisting of particles 
which are allowed to diffuse, compete with one another {2A A), give birth in small areas called 
"oases" {A — > 2A), and die in the "desert" outside the oases 0). This model has previously 

been used to study bacterial populations in the lab and is related to a model of plankton populations 
in the oceans. We first consider the nature of transport between two oases: in the limit of high 
growth rate, this is effectively a first passage process, and we are able to determine the first passage 
time probability density function in the limit of large oasis separation. This result is then used 
along with the theory of hopping conduction in doped semiconductors to estimate the time taken 
by a population to cross a large system. 



I. INTRODUCTION 



Reaction Diffusion Models 



Reaction-difFusion (RD) models have proven to be very 
useful tools for the study of chemical [1], biological [2], 
and ecological [3] systems. RD models typically consist 
of a set of particles which are allowed to diffuse and inter- 
act with one another and their environment in prescribed 
ways. By varying the types of allowed reactions, number 
of types of particles, and reaction rates, one can obtain 
a wide variety of behavior. Much work has been done 
to examine the phase transition between active (popula- 
tion survives as t ^ oo) and absorbing (population dies 
as t — > oo) states [4-8] and to determine the nature of 
propagating fronts [9-12]. 

Typically, RD models are governed by a microscopic 
master equation [13] which describes the probability 
flow into and out of the microstates of the system. 
This master equation is not solvable for all but the 
most simple models, and thus various approximation 
techniques — Langevin equations, for example — are usu- 
ally used. There does exist a systematic expansion of 
the master equation [13], the lowest order of which is 
usually a deterministic differential equation or Fokker- 
Planck equation for the mean concentrations of the con- 
stituent particles. These equations — reaction-diffusion 
equations — are often studied first as a means of char- 
acterizing the qualitative behavior of the model under 
examination; they constitute a mean-field theory for the 
model. 

The effects of quenched disorder in the reaction rates 
on the critical behavior of RD models have been dif- 
ficult to determine. A straightforward renormalization 
group treatment leads to runaway flows [8, 14], but some 
progress has been made using simulations [15-17] and 
real-space RG methods [18, 19]. Disorder effects on RD 
fronts have also been studied, mostly for the case in which 



the disorder is time-dependent ( "annealed" ) and the sys- 
tem admits a front solution in the absence of noise [20-23] 
. However, a few studies have been made of the effects of 
quenched disorder on RD fronts [24] , and some attention 
has been devoted to the interesting case of noise-induced 
fronts [25, 26]. 



B. Our Model: Oases and Deserts 

This work concerns the nature of transport in a partic- 
ular reaction-diffusion system with spatial inhomogeneity 
in the reaction rates. We study a model with a mean-field 
limit defined by the generalized Fisher/KPP equation 



dc{x, t) 
dt 
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D\7^c{x, t)-v-Vc{x, t)+U{x)c{x, t)-qc{x, , 

(1) 

where c[x^t) represents the population density, D is the 
diffusion constant, v is a spatially uniform convection 
velocity (representing the flow of some liquid in which 
the particles exist), U{x) is a spatially inhomogeneous 
growth term fixed in time, and q = Mq'^ is a competition 
term {b is a competition rate and £o is the microscopic 
length scale at which two particles will compete with one 
another). One of the simplest cases to consider is when 
U{x) — —z everywhere except a small patch near the 
origin, where V(x) = y. The region of positive growth 
rate near the origin is called an "oasis," while the rest of 
space is termed the "desert." This model was previously 
studied by Nelson and coworkers [27, 28], and a micro- 
scopic model (the contact process with disorder) with 
this mean-field limit was studied by Joo and Lebowitz 
[29]. Both sets of researchers found a transition in the 
(JJ{x))-\v\ plane between extinct, locahzed, and delocal- 
ized phases in finite systems with periodic boundary con- 
ditions: for high average growth rate and high convection 
velocity, they observed a delocalized phase; for low aver- 
age growth rate and high convection velocity they found 
that the population became extinct; and for low aver- 
age growth rate and low convection velocity they found 
a localized phase. These predictions were tested in a 
laboratory setting using bacteria protected from harm- 
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ful UV light (the "desert") by a mask (the "oasis"); the 
experiments largely confirmed the theoretical predictions 
summarized above [30]. 

In this paper, we will examine the nature of trans- 
port in a system consisting of many identical oases dis- 
tributed randomly at low density in a desert. We term 
this low oasis density regime "hostile" ; the opposite case 
in which oases fill up most of space we call "fertile." 
Because transport between oases in such a system in- 
volves the movement of a low population density, fluctu- 
ations about the mean-field theory (discreteness effects) 
will be important. We will thus be examining a partic- 
ular stochastic process with a mean-field limit given by 



(1). This process is easiest to introduce on a c? = 1 lat- 
tice; the generalization to higher dimensions is trivial. 
Identical particles (labeled A) occupy lattice sites with- 
out occupation number limits and are allowed to undergo 
the following processes: hopping to either side with rate 
w/2 (total hopping rate of w); death {A 0) with rate z 
if in the desert; reproduction {A 2A) with rate y if on 
an oasis; and competition/coagulation {2A A) with 
rate b everywhere. This process is governed by a master 
equation for the joint probability P{{c},t) to have oc- 
cupation numbers {c} = {. . . , c^-i, Cv, Cv+\ . . .} on the 
lattice points u at time t: 



dP{{c},t) 



dt f ( + • • ' ^"^-1 + hcu-l,...,t) + (c,+i + 1)P(. . . , - 1, c,+i + 1, . . . , i) - 2c,P({c}, i)) 

+ ^ [(c, + 1)P(. . . , + 1, . . . , i) - c,P({c}, t)] +Y,y. [(c, - 1)P(. . . , - 1, . . . , t) - c,P({c}, t)] 
+ 6^(c, + l)c,P(...,c, + l,...,t)-c,(c,-l)P({c},t). (2) 

1/ 

I 



Here = on the oases and z in the desert, and 
in the desert and y on the oases. 



II. GROWTH NEAR ONE OASIS 



A. Mean-Field Description 



Let us now present a brief outline of this paper: in 
section II, we will examine the nature of growth near 
a single oasis. Because the mean-field equation for the 
steady-state population density is exactly solvable in one 
dimension, we will be able to identify a length scale de- 
scribing the distance away from the oasis at which fluc- 
tuations about the mean-field theory become important. 
We will also briefly discuss in this section the problem 
of extinction. In section III, we will look at transport 
between two oases. By using the fact that the 2A A 
competition process is unimportant far away from an oa- 
sis where the population is low, we will be able to devise a 
simpler model which captures the transport characteris- 
tics of the full model for large oasis separation. In section 
IV, we will finally tackle the problem of transport in a 
system with many oases. By employing an analogy with 
the problem of hopping conduction in doped semiconduc- 
tors, we will estimate the time taken for a population to 
cross a large system. Finally, we offer a summary of our 
results along with some remarks in section V. Much of 
the material in sections (III) and (IV) has been described 
by us in an earlier publication in less detail [31]. 



We begin with a study of the nature of population 
growth near a single oasis in mean-field theory, starting 
with a ID lattice with a single oasis of width 2a lattice 
points centered at the origin. First, we multiply P({c}, t) 
by Cv in Eq. (2) and sum over configurations to obtain 
an equation for the time evolution of the average particle 
concentration {c^){t): 



(c._i)(t)-2(c.)(t)] 
c,){t)-h{c,{c,-l)){t). (3) 



In order to obtain a "mean-field" description of our sys- 
tem, we replace the term {cv{c^ — 1)) with (cy)^. This 
replacement should work well when the population is 
large — i.e., near the oasis — since we would expect the 
relative fluctuations in particle number to be smaller in 
this case. (There arc, of course, more formal ways of de- 
riving the mean-field equation from the master equation. 
See, for instance, Ref. [13].) With this replacement, we 
can write a mean-field equation for c{u,t) = {Ci,){t): 



dc{u, t) w 
dt ^ 2 



[c{i^ + l,t) + c{iy-l,t)- 2c{u, t)] 

z{u)\c{v,t)-hc{y,tf. (4) 



It is easier to consider the continuum version of this equa- 
tion, which is obtained by introducing a lattice spac- 
ing and redefining c{v,t) c{i',t)£o, b q/f-o, 
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and V x/t^. The diffusion constant D is defined as 
Tliis leads to the d = 1 version of (1), with 
U{x) = {y + z)Q{a— \x\) — z, where 9(-) is the Heaviside 
step function. The length scale £o has an interpretation 
in the continuum as the distance within which particles 
compete with one another. 

There are two things we would like to know: first, what 
does the mean-field concentration c(x, t) look like a,s t 
oo? Second, what is the time scale on which a small 
population grows into a substantial population? Solving 
analytically for c{x, t) for all times is not feasible, but it 
is possible to solve for the steady-state t ^ oo solution 
c{x,t = oo) = Css{x) and thus answer the first question. 
This function is given by 



far away from the oasis. To do so, we drop the nonlinear 
term from the mean- field equation (1) under the assump- 
tion that Css{x) is small far from the oasis. This leads to 
the linear equation 



Css{x) = Css(0)-m+sn^ 



/ q\m- 
6D 



\x\ < a 



\m- 



3^ /fx \ 
Css{x) = -^csch^ (2 -a) + Cj \x\ > a, (5) 



where sn(M, k) is a Jacobi elliptic function, n = y^z/D, 
Css(O) is the steady-state population at the origin, C = 
csch.~ ^ {^/2qE^JcL)J3z) (cssia) is the steady-state pop- 
ulation at the edge of the oasis), and m_|_^_ are defined as 



3c,, (0) - 3y/2q ± v/(3y/2g - c,,(0)) {3y/2q + 3c,,,,(0)) 

The constants Css(O) and Css(a) can be found by match- 
ing the solutions and their derivatives at \x\ = a. This 
leads to a transcendental equation for Css(O): 



m+ sn^ 



s(0) 



6D 



/ 3y-2ge„(0) 



(6) 



Numerically, we have found that an excellent approxi- 
mation to Cas (0) is Css (0) ~ (j/ — yc) I q, where yc is the 
minimum growth rate at which the population does not 
die off as t ^ 00 when q = 0. This cutoff can be found 
by solving (1) with q = (see Appendix B), which leads 
to the following transcendental equation for y^. 



2; cot 




(7) 



At large distances from the oasis (|a;| ^ a), Css{x) 



CooS '''^L where Co 



4-f^Css{a)e'"' (7-1 = l-hcsch(C)). 



In the limit of high growth rate — y 00 with all 
other rates fixed — Css{a) 00 and Cqo 6ze'^°'/q. 
For smaller values of y, Css{a) — and thus c^o — can be 
found by first solving for Css(O) using (6) and then us- 
ing the relation (see Appendix A for derivation) Css{a) = 



3y-2qcss(0) 



3{y+z) ■ 

In higher dimensions, we consider a hyperspherical oa- 
sis of radius a. It is not possible to solve exactly the 

t — > 00 nonlinear mean- field equation for d > I, but it 
is easy to ascertain the asymptotic behavior of Css{x) 







zc. 



(8) 



which is valid far away from the oasis. In two dimensions, 
this is solved by Css{x) ~ CocKq {nr), where r = \x\ and 
Kq is a modified Bessel function of the first kind. In 
three dimensions, Css{x) ~ Cooe~'^^ / nr. Because finding 
an exact solution for the entire space (including r < a) is 
no longer possible for d = 2 or 3, we cannot write down 
an analytic expression for the prefactors Cqo in front of 
these asymptotic functional forms. 

The question of the time scale on which a small pop- 
ulation grows into a substantial population has been ad- 
dressed by Nelson and coworkers [27, 28] . They analyzed 
the eigenvalue spectrum of the linearized [q = 0) version 
of (1) and found that the largest eigenvalue Fq is given 
by [28] 



\ = {y + z)f (^VD/a^{y + z)) - . 



(9) 



where f{x) is a monotonically decreasing function of x 
which goes as 1 — n'^x'^ /A for a; <C 1 and l/x"^ for a; ^ 1 
In the limit of large y, then, Fq y, and the time scale 
on which a small population grows up is 1/y. 



B. Fluctuations and Extinction 

It has been known for some time that fluctuations can 
drive a system to extinction even when mean- field theory 
predicts a stable active state. In the case of a continu- 
ous homogeneous system with the same reactions as our 
system — A 2A with rate y, A ^ with rate z, and 
2A A with rate b — there is an active phase only when 
z — y < Tc, where Vc depends on dimension but is less 
than zero for d = 1,2,3 [8]. Mean- field theory, on the 
other hand, predicts an active phase for y > z; fluctua- 
tions drive the critical growth rate up. The disparity be- 
tween mean-field and stochastic behavior is even greater 
in the case of a d = system: mean-field theory predicts 
a. t ^ 00 steady state which is reached for any nonzero 
initial condition, but solving the master equation leads 
to the conclusion that, for any z > 0, the population 
will eventually become extinct [13]. The mean extinction 
time in this case can be calculated exactly as a function 
of y, z, b, and the starting population no, although the 
resulting expression is cumbersome to work with [13]. 

For the case of a single oasis in an infinite desert, it 
seems clear that the population will become extinct as 
t 00 for d = 1,2,3: the finite oasis cannot compete 
with the infinite desert, regardless of how high the growth 
rate y is. For the problem we will be considering, it is 
important that the oases not die out too early, and thus 
we need to know the dependence of the mean extinc- 
tion time on the various parameters of the problem. The 
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ficld-thcorctic tools used to analyze systems with trans- 
lational invarianee are hard to apply to this case, as are 
the various methods (see Ref. [32] for one such method) 
used to analyze d = systems. Nonetheless, we can try 
to place a lower limit on the extinction time. To do so, 
we will return to the lattice case in one dimension; our 
results will be applicable to the continuum case and to 
other dimensions. 

Consider the case of a perfectly deadly desert, z ^ oo. 
This effectively turns our system into a finite system with 
2a lattice points and absorbing boundaries. The "effec- 
tive" death rate is of the order of w, the hopping rate. 
Now consider & d = system with the same birth and 
competition rates which has a death rate of w, the hop- 
ping rate in our original system. Our d = 1 system 
will certainly live longer than this system, on average: 
the number of events needed to extinguish the popula- 
tion completely is much larger. As mentioned above, the 
mean extinction time for this d = system can be calcu- 
lated explicitly, with the result that Tcxtinct ^ e^^, where 
c is a constant, for large y [13]. This suggests that the 
mean extinction time should rise at least exponentially 
with y in our one oasis problem when y is large. By 
choosing a large y, then, we can ensure that extinction 
will not invalidate our results. Prom here on, we will as- 
sume that the growth rate on the oases is large enough 
that extinction is unlikely on the transport time scales in 
question. 



III. TRANSPORT BETWEEN TWO OASES 

A. Transport as a First Passage Process 

Our eventual goal is to understand the transport of a 
population across a system filled with oases at low den- 
sity. The first step towards such an understanding is 
to determine the nature of transport between two oases. 
Consider two oases of radius a in d dimensions. The cen- 
ter of one oasis is located at the origin, and the center 
of the other oasis is located at position R. At t — 0, the 
first oasis is populated and the second oasis is empty. We 
wish to find the infection time — that is, the time it takes 
for a population to take hold and reach a significant level 
on the second oasis. This time can be roughly broken 
into two parts: Ttransit, the time it takes particles from 
the first oasis to reach the second oasis; and ^gi-owth? 

the 

time it takes the population to rise to a significant level 
once the second oasis has been reached. We will assume 
that the first particle to reach the second oasis will re- 
produce and that its offspring will not die out; in other 
words, we will take TtransU to be the first passage time 
(FPT) of the process. This assumption can be satisfied 
in two ways: the first way is simply to make the growth 
rate y of the oases very high. In this case, it is possible to 
estimate how the survival probability increases as y in- 
creases. Consider again the case of a very deadly desert: 
if the particle diffuses off the oasis, it is certainly dead; 



thus, there is an effective death rate of order D/Iq^. For 
the case of a very small oasis, then, a toy model of the 
oasis is a d = system with death rate of D/io^. For 
this case, it is known that the survival probability goes 
like 1 — D/yto' [13], and thus making y very high assures 
that the population will take hold and survive. A second 
way of satisfying our assumption is to seed the oases with 
a second species of particles, B, which interact with the 
A particles via the reaction A + B ^ 2A aA, a, very high 
rate. 

The time Tgrowth that it takes the initially small pop- 
ulation on the second oasis to grow to a macroscopic size 
should go roughly like 1/y for large y, and so choosing a 
large y should also serve to make Tgrowth ^ ^transit • For 
the remainder of the paper, we will assume that y is large 
enough so that this is the case. Note that by taking y to 
be very high, we have done three things: first, we have 
ensured that a small population which reaches a new oa- 
sis grows into a sizable population and does not die out, 
which allows us to identify the first passage time with 
the transit time; second, we have made the time for this 
growth small compared to the transit time; and finally, as 
mentioned in the previous section, we have ensured that 
extinction will only occur on a time scale much larger 
than the one associated with transit. 

Consider the case where the two oases are close to- 
gether: particles from the first oasis diffuse out in a 
front, its amplitude decaying due to the death term in 
the desert and competition effects. However, so long as 
the second oasis is close enough that the edge of the front 
is almost certain to possess many particles (the number 
will vary from realization to realization of the stochastic 
process), the transit time should simply go as R, the oasis 
separation. However, once R is well above some length 
scale we will call i?iin, this is no longer true: the front 
simply does not exist in most realizations of the system, 
as the number of particles present at this distance from 
the first oasis is quite small for all times. In this regime, 
the second oasis is reached not by a front but by a stray 
particle (or some stray particles) that manages to make 
it through the desert; it is essentially a noise- induced 
growth process. i?iin can thus be roughly defined as the 
distance from the oasis at which the large-time average 
concentration falls to l/^o*^- We have already analyzed 
the mean-field equations for the average concentration 
as t ^ oo, and found that, except in d = 1, there are 
no closed-form solutions. In one dimension, setting the 
mean-field t ^ oo average concentration (5) for large y 
equal to l/£o and solving for R\{^ leads to 




In the limit of large z/h, this simplifies to i?iin ~ a -I- 
\jD/z\n{&z/b), where b is q/^Q. If y is smaller, the rel- 
evant length scale will also be smaller. We believe that 
this length scale should be of the same order of magnitude 
in higher dimensions, and so (10) should also provide a 
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rough estimate of -Run for = 2 and d = 3. 

B. A Simpler Linear Model With a Source 

As we move further from the first oasis, the competi- 
tion process 2A ^ A becomes less and less important, 
especially if b is small compared to the other rates in the 
problem. Due to this fact, it is natural to wonder if ignor- 



ing these interactions altogether might be the first step 
in the creation of a tractable model with the same large 
distance first passage properties as the full model with 
competition. We will now propose such a model, which 
has been discussed by us in an earlier work [31]: consider 
replacing the first oasis with desert, and then placing a 
point source in the middle that produces non-interacting 
particles at some average rate g. The master equation 
for this process on a lattice in d = 1 can be written as 



dP{{n},t) 
dt 



- ( J2(n^-^ + 1)P(. . . , n^_i + 1, - 1, . . . , t) + (n^+i + 1)P(. . . , - 1, n^+i + l,...,t)- 2n^P({n}, t)j 

(11) 



+ zJ2 linu + 1)P(. . . , + 1, . . . , - n,P({n}, t)] + g [P{. . . , no - 1, . . . , i) - P({n}, t)] 

V 

I 



where P({n},t) is the joint probability to have occupa- 
tion numbers {n} = {..., n^-i,n^, n^+i . • •} on the lat- 
tice points V at time t. For an appropriately chosen g, 
the mean flux of particles past the surface at R^^ should 
match that of the model with competitions; beyond that 
point, the model with a source differs from the model 
with competitions only in that it ignores the rare annihi- 
lation interactions between particles. We will show that, 
for an appropriately chosen g, this model which we will 
refer to as the linear model with a source — accurately 
captures the first passage properties of the full nonlinear 
model with competition. 

As with the full nonlinear model with competition 
(hereafter referred to as the nonlinear model), it is useful 
to analyze the mean-field behavior of the linear model 
with a source. The master equation (11) can be multi- 
plied by Ui, and summed over configurations to obtain an 
equation for the time evolution of the average number of 
particles n{v,t): 

' = -[n{y + l,t)+n{v-l,t)-2n{v,t)] 

-zn{u,t) + gS^fi. (12) 

We will study the continuum version of this equation in 
detail in one, two, and three dimensions. Taking the 
continuum limit of (12) (and changing 9| — > for d > 
1) results in: 



dn{x,t) 



DV'^n{x, t) - zn{x, t) + gS'^{x). (13) 



Unlike the mean-field equation for the model with com- 
petitions, this equation can be solved exactly in all di- 
mensions. If we assume an initial condition with no par- 
ticles present, a Laplace transform in time and Fourier 
transform in space leads to: 



n{k, s) = 



9 



s{s + Dk^ + z) ■ 



(14) 



Transforming back into the time domain gives: 



n{k,t) 



Dk'^ + z 



(15) 



We are interested in the long-time, steady-state behavior 

in all dimensions. Letting t ^ oo and transforming in 
space gives the following solutions for nas{x) = n{x,t = 
oo): 



fissir) = 
nss{r) = 



-k\x\ 



gKainr) 



4nD 



AnDr 



ID 
2D 
3D 



(16) 



There is one additional case of interest: the d — 1 
lattice case. The relevant mean-field equation in this 
case is simply (12). After a Laplace transform, we are 

left with a difference equation which can be solved with 



the ansatz n{v + 1, s) 
solution is 

n{y, s) 



-fis) 



n{iy,s) for v > 0. The 



ge 



sw sinh(/(s)) ' 



(17) 



where /(s) = cosh~^(l-|-(s-|-2;)/w). We can immediately 

get the t ^ oo behavior of n{i', t) from this expression by 
multiplying by s and letting s 0, resulting in 



g^ 

nss{i') = n{u,t ^ oo) = — . , (18) 
'u;smh(/) 

where / = /(O). 

The functional forms of the continuum solutions in (16) 

are the same as those of the solutions for the asymp- 
totic (r ^ a) steady-state nonlinear [h ^ 0) equations 
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discussed in Section II A. For a properly chosen cre- 
ation rate g, the mean-field solutions of the two models 
should match at long distances. We will use this method 
of matching mean-field solutions to determine g for the 
purposes of making numerical predictions of first passage 
properties in the nonlinear model. It is important to note 
that g is not a "fit parameter" : its value is completely de- 
termined by the oasis size, the death rate, etc., and is not 
adjusted to fit data generated by the nonlinear model. 

In practice, one can solve the nonlinear steady-state 
mean-field equations numerically, and then find g by 
matching the long-distance behavior to the appropriate 
solution from (16). It is possible, however, to match the 
d = 1 solutions analytically: using the results of Section 
II A together with (16) results in 



(19) 



where 7^ = 1 + yj2qcss{a) /"iz, as before. The con- 
stant Css(o) can be found as described in Section II A. 
As y — > 00, 5 ^ \2\/Dz^ e'^°'/q. For higher dimensions, 
it is necessary to numerically solve the mean-field equa- 
tions for the nonlinear model to accurately calculate g. 



Analytic Predictions from the Linear Model 
with a Source 



With a method in place for determining g from the 

parameters of the nonlinear model, it is now possible to 
use the linear model with a source to make predictions 
about first passage properties of the two oasis system. 
We begin by noting that, since the particles in the linear 
model with a source are non-interacting, the full multi- 
particle FPT PDF /at (a;, t) that is, the probability per 
unit time that the first particle from the first oasis reaches 
the second oasis between t and t + dt — can be written in 
terms of the one-particle FPT PDF .fi{x,t). (Note that 
the vector a; is a stand-in for all the geometric particulars 
of the system. For instance, for a spherical or circular 
oasis, fN{x,t) depends on the distance of the center of 
the oasis from the origin R and the radius a of the oasis. 
These geometrical particulars are not important for our 
present discussion, and so we express /n as a function of 
the generic vector x.) This is accomplished as follows: 
assume the source is at the origin, and that it releases 
N particles per unit time At [33]. Define S{x,t) = 1 — 
Jo dt' fi{x, t') = 1 — Phit(a3, t) to be the probability that 
a particular particle released from the origin at i = 
has not reached the target oasis by time t. If we define 
Pnone{x,t) to be the probability that no particles from 
the source have hit the target oasis by time t, then 



-^none 

{x,t) = JJ [S{x,t)] 



N 



(20) 



Taking the logarithm of this expression gives 

t 

ln[P„one(a;,i)] = Yl 5AHn[5(a;,r)], (21) 

r=0,At,... 

where g = N/At is the creation rate. Taking the limit 
At — > with g fixed and exponentiating both sides leads 
to a closed equation for Pnone{x,t) in terms of S{x,t): 

fnone(a;, t) = cxp (^g dt' In S{x, t')^ . (22) 

Since we are interested in oasis separations large enough 
that a given single particle has a low probability of ever 

reaching the second oasis, S{x,t) is close to 1 even as 
t ^ 00. This allows us to approximate lnS'(a;, t) = ln(l — 
Phit{x, t)) by —Phit{x, t), leading to a simpler expression 

for Pnone{x,t): 
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dt'it-t')Mx,t') . (23) 



The full FPT PDF fN{x,t) is simply -dtPnonc{x,t). 

There is one more useful way to write Pnone: since the 
integral appearing in the exponent in (23) is a convolu- 
tion of t and fi{x,t), its Laplace transform is simply a 
product of the two functions' individual Laplace trans- 
forms. Explicitly: 



P„o„e(a;,i) ^ exp (^-g£ ^ /i(a;, s)/s^ ) , 



(24) 



where C~^[u{s)] is the inverse Laplace transform of u{s) 
and fi{x, s) is the Laplace transform in time of fi{x, t). 
Often it is easier to compute fi{x,s) than fi{x,t), and 
in these cases (24) can be very useful. 

In order to make predictions using (23) or (24), it is 
necessary to compute the one-particle FPT PDF ,fi{x, t). 
We will do this now for the continuum case in all relevant 
dimensions and the lattice case in d = 1. We will start 
with the continuum case. The diffusion equation gov- 
erning the probability distribution pi{x,t) of a particle 
released into the desert from the origin at t = is 



dpi{x,t) 

dt 



= DV'^pi{x,t) - zpi{x,t), 



(25) 



with boundary condition pi (oasis surface, t) = 0. This 
boundary condition is of course not true in the model — 
particles arriving at the oasis will not immediately die — 
but it is used as a device to extract first passage proper- 
ties. By writing pi{x,t) = (f)i{x,t)e^^^ , it is possible to 
eliminate the death term in (25) and arrive at a simple 
diffusion equation for (j)i{x,t). The FPT PDF fi{x,t) 
can be obtained by considering the flux of probability 
into the oasis [34]: 



=o,At,. 



fi{x,t)=D [ 



oasis 
surface 



dA h ■ V0i(a;, t)e 



(26) 
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where dA is an element of the oasis surface and n is a 
unit vector pointing out from the oasis. Since 
is the solution to a simple diffusion equation, D J dA n ■ 
V0i(a;,t) = /f=°(a7,t), the FPT PDF in the case where 
there is no desert. This fact can be combined with (26) 
to arrive at the conclusion 



(27) 



The Laplace-transformed FPT PDF fi{x,s) is thus re- 
lated to the z = function by 



(28) 



These results are convenient due to the fact that, for 
circular or spherical oases, exact solutions exist for 

In one dimension, /f^"(a;,i) = |x|e~^ /y/iirDt^ 
[34]. This means that 



c^o.6 




600 



FIG. 1: Main window: plot showing Pnonc{x,t) in d = 1. 
The lines represent, from left to right, the function for x — 
16, 18, 20, 22, and 24. Inset: a blowup showing the early-time 
behavior of Pnonc- 



y4Dt -zt 



V4^i3/2 



(29) 



when there is a desert present. Plugging this into (23) 
and doing the integration [35] gives 



-fnonc ('^^ ^) 



exp 



4z 

— K|a;| — 



(e''I^IC+ erfc(C+/\/4lI) 
rerfc(r/V4lI))] , (30) 



where C — C i^, t) 



I ±2zt. This function is shown 



in Fig. 1. For large times, Pnonc(a;,^) ^ exp [—ge ''1^1^). 
The j-th moment of fpj{x,t) is given by {T^(x)) — 
j dtPao^c{x, t)P^^; although it is not possible to per- 
form this integral analytically, we can extract its |x| —^ oo 
(large oasis separation) behavior (see Appendix C): 



oK\x\j 



(ID continuum) 



(31) 



In two and three dimensions, it becomes more conve- 
nient to solve for fi{x,s) and then use (24) to obtain 
^rionc- The single-particle FPT PDF is a function of the 
separation of the center of the target oasis from the origin 
(i?) and the radius of the oasis (a), so we will from now 
on write it as fi{R,a,t), where R = \R\. The FPT PDF 
in frequency space in the absence of a desert {z — 0) is 
known for these cases [34]; using (28) gives 



A(i?,a,.) = (^) —— 



K. 



d/2-l 



(32) 



where Kn is the n-th order modified Bessel function of the 
first kind. This equation also holds in d — 1; redefining 
X — R — a and using the definition of K1/2 leads to the 
Laplace transform of (29). 



In d = 2, using (24) and (32) gives 



[R, a, t) ~ exp 



2m 



ds ■ 



'^R 



(33) 

The exponent can be reduced to an analytic function plus 
a real integral (see Appendix C) 



t Ka {kR) 



Y{R,a,t) 



(34) 



RKo (na) Ki (kR) - aKo (kR) Ki (ko) 



ViDi [Ko iKa)Y 



2R'e 



2„-zt r°° 



nD 



du ■ 



where 



Zo ^) 



Jo(f«)2 + yo(»2 



(35) 



The large t behavior of Pnonc is given by Pnonc(^, a, t) ~ 
exp (—g [Kq{k,R)/ KQ{Ka)]t). The moments asymptoti- 
cally approach 



{T^{R,a))=j 



■ \gKo{KR) 



(2D contin.) (36) 



as i? ^ 00. 

The three dimensional case is easy to treat. Since 
K^n{z) ~ Kn{z), looking at (32) immediately shows 
that fi{R,a,s) for d = 3 is identical to the d — 1 
case save for a factor of a/R. Making the replacements 
— > i? — a and g ga/R in (30) gives Pnone{R,a,t); 
making the same replacements gives the t ^ 00 decay 
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PnonciR,a,t) - exp (-5(a/i?)e-''(^-'^)t). The moments 
approach 

/ 7?\ -5 pK,{R-a)j 

(T^(i?,a))=j! - — (3Dcontin.) (37) 

as i? ^ oo. 

The final case we will consider is the d = 1 lattice case. 
Recall that for this case, w is total hopping rate and the 
integer ly denotes the lattice point. The single-particle 
FPT PDF fi{iy,t) is [34] 



t 



(38) 



where is the v-th order modified Bessel function of 
the first kind. It is more convenient to use the frequency 
space function: 



hi'', s) 



' + W + Z+ y/{s + Z){S + Z + 2w) 



(39) 



□ Linear Model with Source 
+ Monte Carlo Data 




200 400 600 800 1000 1200 1400 1600 1800 2000 
Time 

FIG. 2: Binned FPT probabilities for ly = 27 from both the 
linear model with a source (blue boxes) and Monte Carlo sim- 
ulations of the model with interactions (red lines). The width 
of each bin is 50/w, where w is the total hopping rate. The 
error bars on the simulation data represent sampling error. 
Only times up to t = 2000 are shown for the sake of clarity. 



Using this together with (24) gives an expression for 
Pnonc{v,t) (see Appendix C): 



y(z.,t) 



w sinh(/) 



(40) 



e -(«'+^)* r sinie) sin(|z^|6')e^"* ''°<'^'> 
do ■ 



TTW 



[1 



my 



This function decays as Pnonc{'^,t) ~ exp (— ge~-''l''lf) as 
t oo. As in the continuum case, Pnonc cannot be inte- 
grated analytically, but an asymptotic analysis (see Ap- 
pendix C) shows that, as — s- oo, 



(ID lattice) 



(41) 



The agreement between the predictions from the linear 
model with a source and the Monte Carlo simulation re- 
sults from the model with interactions is excellent. The 
linear model with a source correctly predicts the lower 
moments of fM(y,t) for large i/, as shown in Table I. 
A more stringent test of the power of the linear model 
with a source is a comparison of its prediction for the full 
FPT PDF with simulation results. To do this compari- 
son, we integrated fN{i>,t) from (m — l)At to mAt for 
m — 1, 2, 3 ... M to obtain a set of probabilities P(i^, m) 
for hitting the point v for the first time in time bin m. 
We then compared this prediction with simulation re- 
sults. The comparison is shown in Fig. 2 for v = 27; it 
seems clear that the linear model with a source correctly 
predicts the form for /a? {v, t) . 



IV. FROM TWO OASES TO MANY 



D. Simulation Results 

In order to test the predictions of the linear theory with 
a source, we wrote a kinetic Monte Carlo simulation of 
the model with interactions. While it is certainly possi- 
ble to simulate the continuum model in any dimension 
either by doing a discrete-space simulation and choosing 
very small lattice spacings or by using an event-driven 
algorithm [36], we found it more expedient to do a lat- 
tice simulation in c? = 1 and compare with the predictions 
from the lattice version of the linear model with a source. 

Each simulation run began with a population ofy/2b — 
125 particles sitting in the middle site of a 5-lattice-point- 
wide oasis. Once a given end point was reached for the 
first time, the run ended. In order to minimize sampling 
error, 5000 runs were performed. 



A. The Connection with Hopping Conduction 

We have shown that the first passage time statistics 
of the two-oasis model with competition {2A A) are 
adequately captured by a simple solvable model without 
competition when the oasis separation R is large. We 
would like to apply these results to systems with more 
than two oases in order to determine the nature of of 
transport in a large system. 

For concreteness, consider a continuum system in d di- 
mensions {d > 1) comprised of identical oases of radius 
a and growth rate y placed around randomly distributed 
points with number density n in a desert of death rate 
z. We are interested in the low density regime; that is, 
the regime in which the average distance between oases is 
larger than the lengthscale i?iin identified in (10) [37]. We 
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TABLE I: Comparison of predictions from the linear model with a source for the first, second, and third moments with Monte 

Carlo data from the model with interactions. The quoted errors represent a 95% confidence interval. 



Distance 




{T )sim 




{T )sim 




{T )sim 


J/ = 10 

u = 15 
v = 2Q 
iy = 25 
1/ = 30 


12.7781 

33.0596 
102.398 
609.336 
5164.80 


12.4059 ± .092208 

33.3945 ± .288245 
103.966 ± 1.72726 
612.667 ± 15.7994 
5066.48 ± 140.589 


172.986 

1196.65 
14321.6 
6.79632 X lO'^ 
5.26790 X lO'' 


164.97 ± 2.38404 

1223.91 ± 21.8565 
14691.1 ±572.021 
(7.00186 ± .424694) x lO'^ 
(5.13892 ± .321111) x lO'' 


2464.79 

47264.1 
2.69537 X 10** 
1.13187 X 10'^ 
8.05889 X 10" 


2328.16 ± 50.5446 

48850.9 ± 1434.07 
(2.74028 ± .195672) x lO" 
(1.21234 ± .148632) x 10^ 
(7.88352 ± .918500) x lO" 



will allow the oases to overlap, although this shouldn't 
happen too often at the low oasis densities we are consid- 
ering. We will start with one or more oases populated at 
t = and wait for a particular oasis or one of a number 
of oases situated far away to become populated. We will 
call the total time for this to take place Tinfection, the in- 
fection time. Because of the exponential dependence of 
the mean FPT on oasis separation for large oasis separa- 
tions (sec (31), (36), and (37)), the time taken to cross 
the largest oasis separations (or links) on the path should, 
on average, be much greater than the time taken to cross 
the shorter links. The situation is somewhat analogous to 
that of hopping conduction in doped semiconductors [38] : 
the oases in our system play the role of the impurity sites 
in the semiconductor, and the mean transit time between 
oases is akin to the resistance between impurity sites. In 
doped semiconductors, the resistance between impurity 
sites depends exponentially on their separation like e"^, 
where R is the impurity separation and a = 2/a, where a 
is an effective Bohr radius describing the width of the im- 
purity wavefunctions [38] . This is similar to the way the 
mean transit time (and, indeed, all other moments of the 
distribution for large separation) depends exponentially 
on oasis separation in our system. There are a couple 
of significant differences between the two systems: first, 
there is no equivalent in the semiconductor problem of 
the growth time, the time needed for the population on a 
newly inhabited oasis to rise to a significant level; second, 
the resistances between impurity sites are not the aver- 
ages of stochastic variables like the mean transit times, 
but rather definite quantities. The first of these differ- 
ences is insignificant since we have already assumed that 
^growth is much smaller than a typical value of Ttransit for 
oases separated by a large distance. The second differ- 
ence is more important, and some of its implications will 
be discussed in detail later in this paper. 

The problem of determining the resistivity (or conduc- 
tivity) of a doped semiconductor in the hopping regime 
was first tackled satisfactorily using ideas from percola- 
tion theory by Ambegaokar and coworkers [39, 40]. They 
found that the resistivity is dominated by the largest 
links in the network of impurity sites spanning the sys- 
tem. Any links with much larger resistances are effec- 
tively shunted by the smaller resistances, and are not im- 
portant in determining the macroscopic resistivity. The 
size of the largest link i?max can be determined using con- 
tinuum percolation theory, which works in roughly the 



following way: a circle (or sphere) is drawn around each 
impurity site, and the radius of each circle is increased. 
When an impurity site center comes within the circle cen- 
tered around another impurity site, the two are said to 
be linked. When the radii of the circles are increased to 
the point where a cluster of linked sites connects one side 
of the system to another, we have reached the percola- 
tion threshold. The last link to form is clearly the longest 
link, and we call its length i?max- This length varies from 
sample to sample, but has a well-defined limit as the sys- 
tem size goes to infinity [38]: 



Rr, 



Baid) 



nVd 



l/d 



(42) 



where Bc{d) is the dimcnsionally-dependent bonding cri- 
terion, Vd is the volume of a d-dimensional unit hyper- 
sphere, and n is the number density of impurity sites. 
The quantity Bc{d) has an interpretation as the mean 
number of connected neighbors for members of the per- 
colation cluster [38]. 

The network which carries the majority of the current 
in doped semiconductors is called the critical subnetwork, 
and its correlation length is Lq (this is also the length 
scale at which sample-to-sample variations in a-Rmax be- 
come relatively small, of order 1) [38]. Above this length 
scale, the system can be regarded as homogeneous, and 
so the resistivity of a large system of size L ^ is 
roughly equal to the resistivity of a system of size Lq. 
As argued above, this resistivity is largely determined by 
the resistance across the largest link, which is equal to 
6"^"="= /Go (Go has units of conductance and a is equal 
to the inverse). The resistivity p is then given by [38, 40] 



d-2 



Go 



(43) 



with the correlation length Lq given approximately by 

(ai?max)'' 



(44) 



where n is the number density of impurity sites and v 
is a critical exponent equal to 4/3 in = 2 and .9 in 

d = 3. 
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B. Dynamics of Transport in a Macroscopic System 

Now let us return to our problem. Consider a system 

of the size of the correlation length Lq of the snbcritical 
network with one oasis initially infected at one edge of 
the system. In the hopping conduction problem, the goal 
is to find the resistance between the edges of the system; 
in our problem, it is to find the first passage time between 
the starting oasis and either a specific oasis on the oppo- 
site edge or any oasis in a thin layer close to the opposite 
edge. Unlike the hopping conduction problem, our prob- 
lem is dynamic in nature; an additional difference is that, 
as mentioned previously, our first passage times are ran- 
dom variables with a distribution whose mean increases 
exponentially with oasis separation rather than fixed re- 
sistances with an exponential dependence on link size. 
The mean FPT across the system is thus an average of a 
minimum: for a fixed set of oases, each realization of the 
dynamic process yields a path with minimal first passage 
time which may differ from the paths from other real- 
izations. However, there is at least one large link of size 
— ^max which must be crossed in order for the popula- 
tion to reach the opposite edge of the system, and the 
time to cross this link sets the time scale to cross the 
system in the same way the that the resistance of the 
largest link sets the scale of the resistance in the hopping 
conduction problem. Thus, 

(time to cross system of size Lq) — (T(i?max, a)), (45) 

where T(i?inax,a) is given by (36) or (37) depending on 
the dimensionality of the system. 

Now consider a very large system. We wish to find 
the mean infection time (Tinfection(-^^)) — that is, the mean 
time for the population to travel between oases separated 
by some large distance L ^ Lq. This time is roughly 
equal to the mean FPT in the parameter regime in which 
we're interested (that is, the limit of high growth rate 
on the oases). In order to do this, we need to know 
something about the large-scale structure of the cluster 
of oases which will carry the bulk of the particle current. 
Again, looking at the hopping conduction problem is in- 
structive. In that problem, the links-nodes-blobs picture 
[38, 41] suggests that the current-carrying cluster can be 
thought of as a network of nodes separated by a distance 
on the order of Lq connected by one-dimensional links 
and clusters (or blobs) of links. Since the resistance of 
a link depends exponentially on its length, the largest 
one-dimensional links of approximate size i?max largely 
determine the resistance between nodes, and thus the re- 
sistivity of the system, as noted in the previous section. 
(There is some debate as to whether there exists another 
length scale I which, together with Lq, characterizes the 
structure of the current-carrying cluster. See Ref. [42] 
for a discussion of this problem.) 

As a first approximation, let us consider our system as 
consisting of nodes placed on a hypercubic lattice with 
lattice spacing Lq with one large link of size i?max in be- 
tween each node. We ignore the time to cross the shorter 



links and the variations in the oasis configurations from 
one correlation-length-sized chunk to another. The pop- 
ulation starts at one node, and we seek the first passage 
time to some distant node located a distance L away 
along a lattice basis vector (or, equivalently, n = L/Lq 
lattice points away). This is the basic problem of first 
passage percolation (FPP), a field largely studied in the 
mathematical community [43]. One of the basic results 
of FPP is that, as the separation between nodes n — > oo, 
the FPT (r(n)) divided by n goes to a constant /x, con- 
ventionally called the time constant. Thus, the mean 
FPT rises linearly with distance between sites, indicat- 
ing that the proper intensive quantity for our problem is 
the mean FPT divided by oasis separation; in the doped 
semiconductor problem, the proper intensive quantity is 
the resistivity. The value of /x depends on the underlying 
FPT probability distribution, but a general result is that 
M < (^i); where (Ti) is the average time to cross one link 
[43]. (That is, (Ti) is the mean of the distribution from 
which FPTs are picked for each link between nodes.) For 
the case where the times are chosen from an exponen- 
tial distribution, /x ~ -MTi) in two dimensions [44]. In 
general, (Ti) is an upper limit on /i [43]. 

Since we are interested in obtaining a rough estimate 
of the infection time, we will simply use the upper limit 
(Ti) (the mean time to cross one link) as an estimate for 
(Tinfection)/". This givcs us the following: 



(T'infection(-^^)) ^ (T'(-Rmax, «)) 

L ~ To ■ 



(46) 



where (r(i?max5 a)) is again given by (36) in d = 2 and 
(37) in d = 3, and Lq is given by (see (44)) 



(47) 



This result is an order-of-magnitude estimate, but it 
should capture the dependence of (Tinfection) on the rele- 
vant parameters of the system. 

It is probably good to stop at this point and briefly 
recall the approximations that we have made to obtain 
the result in (46): first, we have ignored the growth time 
on the grounds that it is small compared to the transit 
time between oases; second, we have simplified the pic- 
ture of transport on the scale of Lq, replacing the mess of 
oases with a single link of size i?max; and third, we have 
used an upper limit on the time constant rather than the 
time constant itself. It should be noted that the first and 
second approximations tend to lead to underestimating 
(^Infection), while the third tends to lead to overestimat- 
ing it. 



C. Comparison with Simulations 

In order to confirm the predictions of the preceding 

section, we wrote a program capable of simulating a very 
large system in two dimensions. To make the simulation 
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of such a large system tractable, we made some impor- 
tant simplifications which must be explained. The first of 
these is the most important: rather than simulating the 
motion of individual particles, we simply assigned first 
passage times between oases. This allowed us to go to 
system sizes many orders of magnitude larger than we 
could have achieved via a full kinetic Monte Carlo simu- 
lation involving every particle. 

The second simplification involves the nature of the 
FPT PDF used to generate the passage times between 
oases. The linear theory with a source produces an ana- 
lytical expression for this FPT PDF (see Eqs. (34) and 
(35)), but this is unwieldy and computationally expen- 
sive to calculate. However, for large R, the moments of 
this FPT PDF in d = 2 approach those of an exponen- 
tial distribution with parameter gKo{KR)/ Ko{Ka) (see 
36), where k = sJzpD. Since it is the large- i? sepa- 
rations which will largely determine the infection time, 
we simply replaced the complicated FPT PDF between 
oases with this exponential distribution; the errors in- 
troduced by this simplification are serious only for small 
oasis separations, and these do not contribute much to 
the infection time. 

The remaining simplifications are minor: we treated 
all the oases as points; we ignored the growth time, just 
as we have done in the analytical work presented in the 
preceding sections; and finally, we ignored the effects of 
neighboring oases on the first-passage time statistics be- 
tween two oases. This final simplification again intro- 
duces errors mostly in areas of high oasis density where 
oasis separations are small. The bottlenecks of our par- 
ticle current- carrying cluster occur where there are two 
oases separated by a large region of desert, and in these 
areas the FPT statistics should be very close to those 
derived in the case of two oases in an infinite desert. 

Before presenting our simulation results, we must first 
provide some details of the way time was scaled in our 
simulations. With the simplifications we have made, the 
FPT PDF f]\[{R,a,t) between a pair of oases of radii a 
separated by a distance R can be obtained from (34). It 
is given by: 



fN{R,a,t) = — — — - cxp 

Kq (Ka) 



9K0 (kR) 

Kq (kg) 



(48) 



If we define the variable r = giKaiKR^i^y^)/ Ko{Ka) — 
effectively measuring time in units of the mean time to 
cross a link of size -Rmax — we can absorb the dependence 
of the FPT PDF on a and g into r. The FPT PDF then 
becomes 



fN{R,T) 



Ko{kR) 
ii'o(Kiimax) 



exp 



Ko{kR) 

-K'o(K-Rmax) 



(49) 



Our simulation measured time in units of r, so that there 
was no need to input information about g or the oasis size 
a. 

There is one further approximation that we made in 
our simulations simply for the sake of convenience: we 



used the large- argument asymptotic form for Ko{x) of 
^7r/2a;e-^, making the FPT PDF 



fNiR,T) 



R 



exp 



R 



max n{R-R„ 



R 



(50) 

Like some of the other simplifications and approxima- 
tions we made in the simulations, this approximation is 
not good for small oasis separations, but the errors intro- 
duced arc ultimately unimportant given the contribution 
of the small oasis jumps to the transit time. 

If our theory is correct, the mean time to cross one 
block of size Lo in these units (in units of r) should be 
of order 1, and the mean infection time should be 



(Tinfection(-^)) — ( "J^ 



(51) 



If K and -Rmax are adjusted in such a way so that their 
product remains constant, then this amounts to a triv- 
ial rescaling of space, and rinfoction should simply vary 
as 1/-Rmax- This is already captured through the depen- 
dence of Tinfection(-t') ou Lq, and SO we can rewrite (51) 
as 



('rinfectioii(i)) = ) F{KRn 



(52) 



where F{KRmax) is some function of order unity. We thus 
expect that a graph of (rinfoction) versus L/Lq for large 
L should be a straight line with slope of order 1. 

For each simulation run, k and the oasis density n were 
input, -Rmax and Lq were calculated from (42) and (47), 
respectively, and a starting oasis was chosen near the 
center of the system. The simulation then proceeded one 
infection event at a time, with infection times between 
oases generated using the distribution given in (50). In 
order to speed up the simulation, we set a maximum dis- 
tance -Rcut beyond which oases were effectively discon- 
nected. This allowed us to generate new oases "on-the- 
fiy" as the simulation proceeded; together with our prac- 
tice of throwing away information about an oasis once it 
was reached, this allowed us to only keep a small subset 
of oases in memory at any one time, thus allowing for the 
simulation of very large systems. The value of -Rcut was 
chosen so as to make the probability of a missed event — 
that is, a jump event of size larger than -Rcut occurring 
over the course of the simulation — very small (< 10~^). 

In early simulation runs, we found that our starting 
oasis would sometimes be isolated from the rest of the 
cluster, leading to larger-than-expected infection times 
with a large contribution from the time for the popula- 
tion to make the first jump. In the limit as -L ^ 00 — 
the large-distance limit we're interested in — this contri- 
bution to the infection time, which docs not grow with L, 
should become negligible, but for finite values of L it can 
be important. In order to eliminate this effect from our 
simulations without going to system sizes too large to be 
simulated in a reasonable amount of time, we allowed the 
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FIG. 3: First passage times across a large system shown for 
seven different combinations of k and 7?max. Error bars are 
not sfiown since they are, in most cases, smaller than the 
symbol size. The lines represent best-fit lines for each k, i?max- 
The two lines with fi;J?max ~ 12.0 lie nearly on top of one 
another, as one would expect; we have omitted every other 
data point for each of these runs for clarity. Note that the 
value of the slope (which is equal to F(Kiimax)) increases as 
tiimax increases. 



population to "find" the cluster: we restarted the sim- 
ulation once an oasis at least 2i?,iiax from the starting 
oasis had been hit with the newly hit oasis as the new 
starting oasis. The choice of 2i?max was admittedly ar- 
bitrary, but it did serve to eliminate the undesired effect 
from our simulations. 

Once the population was restarted, the simulation con- 
tinued one oasis infection event at a time. When an oasis 
within a small distance 5 <C of one of a set of con- 
centric rings centered at the starting oasis was hit, the 
time and distance from the starting oasis were recorded; 
once all oases in some final ring were infected, the simu- 
lation ended. The results of the simulation are shown in 
Fig. 3. The data confirms our picture of transport: the 
slopes of the best-fit lines through the data are indeed 
of order 1, suggesting that i?max is the correct length 
scale of the largest jumps the population must make on 
its way through the system and that Lq is the correct 
length scale for the distance between these large jumps 
(of course, the population left behind the front edge will 
eventually make larger jumps to infect isolated oases, but 
this is unimportant in trying to determine the infection 
time) . Note that there are some "missing" points on the 
two lines with the highest /vi?max- This is due to the pres- 
ence of oases inside those rings which were not hit before 
the simulation time ended. As Ki?max is increased, such 
outlying oases take longer to hit, but since their "extra" 
contribution to the mean transit time does not scale with 
L, they do not affect our L ^ oo results. 

The slope for each line is equal to the scaling func- 
tion -F(Ki?max) for those values of k and iimax! note that 
i^(Ki?,nax) appears to increase for increasing values of 



Ki?inax- This is likely due to that fact that, as Ki?max in- 
creases, the correlation length Lq increases, and thus the 
number of smaller oasis separations between the large oa- 
sis separations increases as well. We do not understand 
this phenomenon completely, but it is seems a good can- 
didate for further study; however, as the slopes are all of 
order 1, an understanding of this phenomenon is hardly 
essential for making our present argument. 



CONCLUSIONS, REMARKS, AND FUTURE 
WORK 



In this paper, we have examined transport in a 
reaction-diffusion system with disorder in the reaction 
rates. Such systems have been used in the past to 
model bacterial population dynamics and the movement 
of plankton in the oceans. Our model consists of parti- 
cles which are allowed to diffuse with diffusion constant 
D and compete for resources {2A A) everywhere with 
rate 6, but which can only give birth (A 2A) on small 
patches called oases at rate y and which die {A ^ 0) ev- 
erywhere else at rate z. We have considered the limit in 
which the growth rate on the oases is very high and the 
oasis density is very low; in this limit, the time needed for 
a small population to grow on an oasis is much smaller 
than the typical time needed to jump from oasis to oasis, 
and thus transport can be thought of as a first passage 
process. Because the population density traveling from 
one oasis to another is small, it is necessary to consider 
discreteness effects. In order to determine the first pas- 
sage time probability density function (FPT PDF) be- 
tween two oases, we have employed a simplified model 
in which competition is ignored and the initially infected 
oasis is replaced by a particle source. Simulations sug- 
gest that this model correctly predicts the FPT PDF for 
large oasis separations. 

We have used an analogy with the theory of hopping 
conduction to argue that the largest oasis separations 
in the particle current-carrying cluster largely determine 
the time taken for a population to travel to a given tar- 
get. The scale of these separations can be found using 
continuum percolation theory, as in the hopping conduc- 
tion problem. There is a significant difference between 
the two problems: ours is dynamic, while the hopping 
conduction problem is not. However, the use of results 
from first passage percolation theory suggest that the 
time scale for transit should still be determined by the 
largest oasis separations in the relevant particle current- 
carrying cluster. 

There are certainly many future areas of study related 
to our work. First off, there is the obvious question 
of what happens when the oases are not identical, but 
instead have their sizes and growth rates picked from 
some distribution. One might hope that the theory of 
variable-range hopping [38] would be useful in this case, 
though it remains to be seen whether the dynamic nature 
of the problem would make a fruitful mapping possible. 
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There is also the problem of determining the nature of 
the front that moves through a system lik(^ tlic^ one stud- 
ied in this paper. The veloeity of such a front should 
be given roughly by Lo/(I"(-Rmax, a)), but its shape is 
an open question. Finally, there is also the more gen- 
eral problem of RD wavefronts in media with quenched 
disorder, which is a challenge for future studies. 
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APPENDIX A: SOLUTION OF THE 
STEADY-STATE MEAN-FIELD EQUATION IN 



The equation we need to solve is 

= DcJ^{x) + U{x)css{x) - bcss{xf, 



(Al) 



where U{x) = {y + z)Q{a—\x\)~ z and the primes denote 
differentiation with respect to x. To solve this, we can 

find solutions in the oasis < a) and desert > a) 
and then match at the boundaries. In the desert, the 
relevant equation is = DcJ^{x) — zcssix) — bcss{x)'^- 
We define u(c) = u{css{x)) = Cg/(a;), which leads to the 
first order equation 



du z b 2 

= u— —c — —c , 

dc D D ' 



(A2) 



where we have written Css{x) as c for simplicity. This 
equation can be integrated to give 



u{c) = 



dc{x) 
dx 




(A3) 



which can in turn be integrated to obtain the function 
Cssix) quoted on the second line of (5). A very simi- 
lar procedure can be done for the area inside the oasis, 
leading to 



u{c) 




C(X)3-C(0)3-|(C(X)2 



c(0)2 



(A4) 

which can also be integrated, leading to the function 
quoted on the first line of (5). Since the derivatives must 



match at the boundary, we can set (A3) and (A4) equal 
at |a;| = a and obtain the following relation: 



Css{a) = Css(0)4 



/3^-26c,,(0) 
3(2/ + 2) 



(A5) 



APPENDIX B: DERIVATION OF THE 
FORMULA FOR j/c 

The cutoff vahic of the growth rate y below which a 
population placed on an oasis will die out as t ^ oo can 
be estimated using the mean- field equation (1) with b = 
0. For values of y greater than the cutoff, the population 
will continue to increase without limit as i — > oo; for 
y < yc, the population will eventually die out. At yc, 
there will be a steady-state solution. Hence, one way 
of finding the cutoff is to try to match solutions to the 
steady-state equation for \x\ < a and \x\ > a at |a;| = a. 
Only along a certain line in parameter space will this be 
possible. 

In one dimension, the steady-state mean-field equation 
with 6 = is solved by c(0) cos{-^/y/Dx) for \x\ < a and 
c(a)e~'*(l^l^°' for \x\ > a. Matching the functions and 
derivatives at I a; I = a leads to: 



Vc 



^;cot 




(ID) 



(Bl) 



which is precisely (7). In two dimensions, a similar cal- 
culation leads to 



Vc 



Jo{^a)K, (Kg) 



Ji{^a)KoiKa) 
while in three dimensions we have obtained 



(2D) (B2) 



yc = z tan 




(3D) (B3) 



These equations can be solved numerically to determine 
yc- A plot of yc as a function of z in one, two, and three 
dimensions, with all other parameters fixed, is shown in 
Fig. (4). 



APPENDIX C: ASYMPTOTIC ANALYSIS OF 
THE MOMENTS OF /jv(x,i) 

In this appendix, we derive the results for the asymp- 
totic moments of fN{x,t) quoted in Eqs. 31, 36, 37, and 
41. We start with the continuum case. In any dimen- 
sion, Pnone(i?, a, f) ~ exp[— ^^(i?, a, f)] , where Y[R,a,t) 
is given by 



Y{R,a,t) 
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FIG. 4: Cutoff growth rate yc as a function of death rate z 
with a — 3.0, D — 0.5. Here aoi is the first zero of Jq. Note 
that in one and two dimensions, an arbitrarily small growth 
rate with z — will allow a stable population to take hold; in 
three dimensions, ydz = 0) = ■k^D/Ao? . 




Rep 



FIG. 5; Schematic of the contour integral which must be done 
to find Y[R, a, t). The dashed lines represent contributions to 
the integral that vanish as they are moved further from the 
origin. 



Although in 0? = 1 and d = 3 this Laplace transform 
machinery is unnecessary — we can simply perform the 
integral over time appearing in Eq. 23 — it is easier 
to determine the asymptotic behavior of the moments 
of fN{R,a,t) in all dimensions by using these tools. 
Changing variables to p = s + z leads to Y{R,a,t) = 



[{a/Rr/ 



2-1^ 



7(27rz)]Qi(i?,a,t), where 



K. 



d/2-l 



Qi{R,a,t) = / dp- — 



(C2) 



This integral can be evaluated using contour integral 
techniques. There is one second-order pole at p — z and 
a branch cut which we will take to lie on the real p axis 
from p = to p — —oo. Our countour will be taken to 
enclose the pole at p = z, and consists of three parts: 
Qi, the value of which we wish to find; and Q2 and Q3, 
whose values must add with that of Qi to equal 27riS, 



where S is the residue at p = z. The space is shown 
schematically in Fig. 5. Using the residue theorem and 
changing integration variables to m = —p gives: 



2m- 



2m 



e^'^'aK^ (kR) K^+i {kg) 



(C3) 



du 



Mfj,{R,a,u) 



(u + z)2 {-^V^ a) 



where we have used /i = d/2 — 1 and Mf^(R, a,u) 



2zlm [K^{i^R)K^ 



We 



that 



Y{R,a,t) thus has the form Cit — C2 + C^h{t), where 
the Cn are constants in time and h[t) is given by some 
complicated integral. Since y(i?,a,0) = 0, we can let 
C3 — C2 and /i(0) = 1. It should be clear that h{oo) = 0, 
and that h{t) < 1 for all t. This is enough to prove 
the asymptotic results for the moments of f]\[{R^ a,t) 
quoted in Section IIIC. These moments are given by 
{T^{R,a)) = j dtPnonc{R,a,t)P~^; plugging in the 
form for Y{R, a, t) gives: 

/"OO 

{T\R,a))=] dte~slCit-C^(i-h(t))]^,-i ^^^^ 



The constant C2 go to as i? ^ 00, so one can Taylor 
expand exp[gC2{t — h{t))] and arrive at 

/>oo 

{T\R,a))^] / die-s<^i*i^"-Ml + 5^2(1 + 



(C5) 

Keeping only the lowest order term, we get {T^{R,a)) ~ 
j!(.gCi) as i? — > cx). Looking at (C3), we see that 
Ci'= {a/Rf/^-^Kd/2-l{^^R)/Kd/2-l{^^CL). We are now 
ready to plug in the functional forms for and 
arrive at the final asymptotic expressions for {T^ (R, a)): 



{T^iR,a)) - jl 



ID 



9' 

gKoiKR) 



(C6) 



2D 



(T^(i?,«)) = .!(-) ^— 



3D, 



where \x\ — R ~ a (the distance from the origin to the 
edge of the oasis nearest the origin) . 



APPENDIX D: CONVECTION EFFECTS ON 
FIRST PASSAGE PROPERTIES 

We wish to determine the effects of a small convec- 
tion velocity on the first passage properties of a system. 
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Physically, such a convection velocity might represent the 
effects of a moving liquid medium in which the particles 
exist. We start with a two-oasis system and use the linear 
model with a source to make analytical predictions pos- 
sible. To begin with, we replace the initially populated 
oasis with a source located at R and center our coordi- 
nate system in the middle of the target oasis of radius 
a. The convection velocity v is taken to be constant in 
space. In order to solve for Pnono(-Rj we must find 
fi{R, a, t), the single-particle FPT PDF. This is done by 
solving for pi{x, t), the probability density function of a 
single particle released from the source at i? at t = 0, 
and then finding the probability flux into the oasis. 
The diffusion equation governing pi (a;, t), , is 

dpi{x,t) _ 2 



dt 



DV^Pi{x,t)-zpi{x,t)~v-ypi{x,t). (Dl) 



It is essential to simplify this equation before proceeding 
with a Laplace transform. As with (25), we can define 
a new function (j)i{x,t) = pi{x,t)e^* and eliminate the 
—zpi{x,t) term from the equation. We can further de- 
fine the function xi(£C,f) via (j>i{x,t) = e"'^/'^^xi{x,t), 
leading to 

= DV\,{x,t) -^Xi{x,t). (D2) 

The last term on the right can be handled by defining 
Xi{x, t) = ipi{x, t)e~'" */'*^, leading to a simple diffusion 
equation for ipi. 

The flux into the oasis can be used, as before, to find 
fi{R,a,t): 

fi{R,a,t) = Da^-^ j dfl drPi{x,t) (D3) 



where dfl is a differential element of angle in 2D, and of 
solid angle in 3D. All that must be done is to find ')pi{x, t). 
This function is the solution to a simple differential equa- 
tion with initial condition '4)i{x, 0) = e~'"'^/'^^ 5'^{x — R), 
and is thus equal to e~'"'^^'^^(j)i{x,t), where (t)i{x,t) is 
the solution to the simple diffusion equation in the ab- 
sence of convection. Thus, 



(D4) 

We are interested in the case where R ^ a, and so a 
decent approximation of /i {R, a, t) is given by 



MR,a,t) ^ e-''^*/^^e-«/^^/r°(i?,a,i). (D5) 



Note that in the above equation, we have reversed the 
sign of R since it is more natural to take the source 
as the origin rather than the center of the target oa- 
sis. This result can be used to determine the moments of 
fN{R, a, t). By making the replacements z ^ z + v'^/AD 
and g gQv R/2D ^j^g expressions for the moments of 
fwiR, a, t), we arrive at the following expression, valid in 
any dimension: 



(T^(i?,a)).=i!(f) ^ 



where = ^/z/D + v'^/AD'^. 
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